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We study Bianchi cosmologies in the ghost-free bigravity theory assuming both metrics to be 
homogeneous and anisotropic, of the Bianchi class A, which includes types I,II,VIo,VIIo,VIII, and IX. 
We assume the universe to contain a radiation and a non-relativistic matter, with the cosmological 
term mimicked by the graviton mass. We find that, for generic initial values leading to a late-time 
self-acceleration, the universe approaches a state with non-vanishing anisotropies. The anisotropy 
contribution to the total energy density decreases much slower than in General Relativity and shows 
the same falloff rate as the energy of a non-relativistic matter. The solutions show a singularity 
in the past, and in the Bianchi IX case the singularity is approached via a sequence of Kasner-like 
steps, which is characteristic for a chaotic behavior. 
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I. INTRODUCTION 

The recent discovery of the ghost-free massive grav- 
ity theory [1] and its bigravity generalization [2] has re- 
vived the old idea that gravitons can have a small mass 
[3] . Theories with massive gravitons were for a long time 
considered as pathological, mainly because they exhibit 
the Boulware-Deser ghost - an unphysical negative norm 
state in the spectrum [4]. However, there are serious rea- 
sons to believe [5] that the theories of [1],[2] are free of 
this pathology, hence they can be considered as healthy 
physical models. 

Theories with massive gravitons can be used in order to 
explain the observed acceleration of our universe [6] . This 
effect can be accounted for by introducing a cosmological 
term to the Einstein equations, however, this poses the 
problem of explaining the origin and value of this term. 
An alternative possibility is to consider modifications of 
General Relativity, and theories with massive gravitons 
are natural candidates for this, since the graviton mass 
can effectively manifest itself as a small cosmological term 
[7] . This justifies interest towards studying cosmological 
solutions with massive gravitons. 

The known cosmologies with massive gravitons can 
be divided into two types. The first type is provided 
by solutions described by two metrics which are not si- 
multaneously diagonal. For these solutions the graviton 
mass gives rise just to a constant term in the Einstein 
equations and to nothing else, at least for the homo- 
geneous and isotropic backgrounds. Such solutions had 
been first obtained without matter [8], while later the 
special [9], [10], [11] and also general ]12] solutions includ- 
ing a matter source were found. For these solutions the 
mater dominates at early times, when the universe is 
small, while later the mass density decreases and the ef- 
fective cosmological term becomes dominant, leading to 



* maeda@waseda.jp 

t voIkov@lmpt.univ-tours.fr 



a self-acceleration. Such solutions exist for all (open, 
closed and fiat) Friedmann-Lematre-Robertson- Walker 
(FLRW) types, both in massive gravity and bigravity. 
However, perturbations around these backgrounds are 
expected to be inhomogeneous - due to the non-diagonal 
metric components, hence solutions of this type are some- 
times called 'inhomogeneous' ]10]. 

For the second type solutions both metrics are diago- 
nal and of the FLRW type. In this case the effect of the 
graviton mass can be more complex and reduces to that 
of a cosmological term only at late times. Such solutions 
have a somewhat narrow existence range. For example, 
in massive gravity they exist only for the spatially open 
FLRW type [13[, and although in bigravity they exist for 
all spatial types [14[,[15[, they do not admit the limit 
where one of the two metrics becomes fiat. The bigrav- 
ity solutions exhibit rather complex features and show 
several branches. There are physical solutions for which 
the matter dominates at early times, while the graviton 
mass becomes essential later. In addition, there are also 
exotic solutions for which the graviton mass contribution 
is dominant at all times. 

Up to now, cosmologies with massive gravitons have 
been studied mainly in the FLRW limit. It may be im- 
portant to extend the analysis to more general space- 
times. For example, one may wonder whether the FLRW 
universe is stable against anisotropic or inhomogeneous 
perturbations [16[,[17[. Another motivation is the so- 
called cosmic no-hair conjecture. In General Relativity 
(GR) with a cosmological term the de Sitter spacetime 
is expected to be an attractor [18[. In massive gravity or 
bigravity the graviton mass gives rise to a cosmological 
term at late times, and so the acceleration of the uni- 
verse can be explained by the de Sitter expansion. How- 
ever, since the graviton mass is not exactly a cosmologi- 
cal term, one may wonder whether the de Sitter space is 
still an attractor for generic initial values. An interesting 
byproduct in such studies can be an observational relic - 
if small anisotropy remains in the accelerating phase, it 
could be observed [16[. 
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In what follows, we undertake a systematic analysis of 
anisotropic cosmologies in the ghost-free bigravity [2] . as- 
suming both metrics to be simultaneously diagonal and of 
the same Bianchi type within the Bianchi class A, which 
includes types I,II,VIo,VIIo,VIII, and IX. 

As a starting point, we find exact solutions for the 
Bianchi I type for which the two metrics have iden- 
tically the same anisotropics, which however requires 
to fine-tune the two matter sources. Then we attack 
the problem numerically and discover that, even in the 
presence of a matter and for all other Bianchi types 
of class A, the equal anisotropy configurations play the 
role of late time attractors. Specifically, starting form 
arbitrary initial data, the anisotropics approach equal 
and constant values. For Bianchi I solutions constant 
anisotropics can be scaled away, but not for other Bianchi 
types, so that generic homogeneous spacetimes run into 
anisotropic states. 

We find that the shears approach zero exponentially 
fast as the universe approaches the de Sitter phase. How- 
ever, since the anisotropics do not approach zero but os- 
cillate around constant values, the shear contribution to 
the total energy density decreases only as an inverse cube 
of the size of the universe. This is the same falloff rate 
as for a non-relativistic matter, whereas in GR shears 
decrease as the inverse sixth power of the size of the 
universe. Therefore, the anisotropy effect could be ob- 
servable, and so it would be interesting to compare our 
predictions with observations. Since the anisotropy con- 
tribution shows the same falloff rate as a cold dark mat- 
ter, it is tempting to think that the latter could in fact 
be the effect of the anisotropics, although it is unclear 
if this interpretation can also explain the dark matter 
clustering. 

The rest of the paper is organized as follows. In the 
next two sections we describe the ghost-free bigravity 
[2], the Bianchi cosmologies, and derive the field equa- 
tions in the case where the two metrics are simultane- 
ously diagonal and of the same Bianchi type of class A. 
In Section IV we study exact solutions of these equa- 
tions for the Bianchi I type. First, these are solutions 
with proportional metrics, they are the same as in GR. 
Next, these are solutions with identical anisotropics for 
the two metrics, which can be of the FLRW type [14], [15] 
if anisotropics vanish. In Section V we present the gen- 
eralization to the other Bianchi types of class A. We de- 
scribe our procedure for the numerical integration - the 
dynamical system formulation of the problem and the im- 
plementation of the initial values constraints, after which 
we present our numerical results. Section VI contains 
concluding remarks, while in the Appendix we describe 
the Bianchi I solutions in GR and also give the details of 
the dynamical system formulation. 

We set h = c = 1 and use the sign conventions of 
Misner-Thorne- Wheeler. The space and time scale is cho- 
sen to be the inverse graviton mass. 



II. THE GHOST-FREE BIGRAVITY 

The theory is defined on a four-dimensional spacetime 
manifold equipped with two metrics, g^^ and f^^. The 
kinetic term of each metric is chosen to be of the stan- 
dard Einstein-Hilbert form, while the interaction between 
them is parametrized by a scalar function of the tensor 

7t - (2.1) 

where g'^" is the inverse of gp,y and the square root is 
understood in the sense that 

{j'r, ^ 7^,7". = g^"fa.. (2.2) 

We also assume a g-matter and an f-matter interacting, 
respectively, only with g^i, and with f^,^. The action is 

^[g, f , matter] 

J ^>^f J 

-"^ Jd'xV^'^igJ] (2.3) 

+ 4'°' [g, g-matter] + S^^^ [f , f-matter] , 

where R and TZ are the Ricci scalars for g^^ and f^i,, 
respectively, = 8ttG and — 8ttQ are the corre- 
sponding gravitational couplings, while and 
m is the graviton mass. The interaction between the two 
metrics is given by 

4 

^ = Y,bk^kh), (2.4) 

n=0 

where bk are parameters, while '^kij) are defined by the 
following relations^ 

'%(7) = -^W-e^''''" = l> (2-5) 
^^1(7) - - | ep.p.e"'""^7''. = E ^-4 - [7], 

A 

A<B 

'^3(7)--^W-e"''""7^a7>^^= E ^^^B^c 

A<B<C 

= ^([7p-3[7][7']+2[7^]), 
^^4(7) = - ^ ep.p.e"^^S^o7''/j7^7''5 = A0A1A2A3 

= Ub]' - 6[7]'[7'] + 8[7][7'] + - 6(7']) ■ 



Notice that e^^^^ = -€0123 = 1. 
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Here {A = 0,1, 2, 3) are eigenvalues of 7'^^. and, using 
the hat to denote matrices, we have defined 

[7] = tr(7) = b"] = tr(7') = il'T, ■ (2-6) 

It is also useful to express the interaction in terms of 



(2.7) 



k=0 



where '^^(/C) are defined by the same expressions as in 
(2.5), up to the replacement — > fj,A = 1 — A^ and 7 — 
/C, where /i^i's are eigenvalues oiK^^,. The parameters Cfc 
are related to bk as 



Co = feo + 46i + 6&2 + 463 + b^, 
ci = -(61 + 362 + 363 + 64), 
C2 = 62 + 263 + 64, 
C3 = -(63 + 64), 
C4 = 64 , 



(2.8) 



with the inverse expressions 



^0 = Co + 4ci + 6c2 + 4c3 + C4, 

61 -(ci + 3c2 + 3c3 + C4), 

62 C2 + 2C3 + C4, 

63 = -(C3 + C4), 

64 = C4. 



(2.9) 



In the weak field limit, when g^i, « f^jy « 77^,^, the tensor 
/C tends to zero while '^k ~ IC^ so that Eq.(2.7) gives the 
interaction in terms of powers of deviation from the flat 
space. We require the flat space to be a solution of the 
field equations, which is only possible if 



Co 



ci =0, 



(2.10) 



while the quadratic part of the interaction should repro- 
duce the Fierz-Pauli term, 



which fixes the normalization 



C2 = -1. 



(2.11) 



(2.12) 



The above choice implies that the bare cosmological con- 
stant is zero, so that a non-zero cosmological term can 
only be of a dynamical origin, due to the graviton mass 
contribution. Terms proportional to C3 andc4 can be 
kept, so that 



60 = 4c3 + C4 — 6, 61 = 3 — 3c3 — C4, 

62 = 2C3 + C4-1, 63 = -(c3 + C4), 

64 = C4. 



(2.13) 



A. Field equations 

Assuming the spacetime coordinates x^^ to be dimen- 
sionless, the metrics g^i,, f^^ have the dimension of length 
squared. To pass to dimensionless quantities, we make a 
conformal rescaling 



Varying the action then gives the field equations 

1 



-.2 — fJ-i^ 



'my I T J— I 



(2.14) 



(2.15) 



(2.16) 



where G^^ and Q^^, are the Einstein tensors for g^^ and 
f^ij. The graviton energy-momentum tensors are 



9 1 2- 

2 \ §gfJ,U 



(2.17) 



In order to perform the variations here, one uses the re- 
lations 



77 77 

2 3M«(7"n- 2 3.o(7X 



71 U 

gf,^ = - 2 UciYT. = - 2 /-(7")";. , (2.18) 

which can be obtained by varying the definition 7^^ = 
Vg^'^faiy and using the properties of the trace. Intro- 
ducing the angle 77 such that 



K cos 77, 



if = Ksmr], 



one obtains dimensionless quantities 

cos2r7(rt-^,5^;). 



^2^[7]a._^ = - sin^ ?7 ■ 



'9 



'f 



(2.19) 



(2.20) 



(2.21) 



where ^ .g^"T„^ and = /^"Ta,. while 

rt - {61 ^0 + 62 ^1 + 63 '^2 + 64 ^^317" 

-{62^0 + 63^1+64'^2}(7')''. 

+ {63^o + 64^i}(7')''. 

- 64 ^0 (7')^. (2.22) 
with = ^i:(7). The matter energy- momentum ten- 
sors r'™|/ and T'™], are dimensionful, we assume them to 
be of perfect fiuid type. However, they enter the equa- 
tions only via dimensionless combinations'^ 

4 = tHm^ ^ (p + Pg)u^'u, + Pg Si; 

^2 

^ 7-[m]M_^ = 7-[m]p^ ^ ^ Py)W'^Z^, + Pfdi;. (2.23) 



^ We include the gravitational constants and and the gravi- 
ton mass m in the definition of the energy densities and pressures. 
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As a result, from now on we shall be considering the 
field equations (2. 15), (2. 16) expressed entirely in terms 
of dimensionless quantities. 

The diffeomorphism invariance of the matter terms in 
the action implies the conservation conditions 



(9) 



(/) 



'Hp 



0, 



(2.24) 



where V and V are covariant derivatives with respect 
to and f^^. The Bianchi identities for (2.15) then 
imply that 



(a 



(2.25) 



Similarly, the Bianchi identities for equations (2.16) im- 



ply that Vu / 



Mm 



0, but in fact this condition is 



not independent and follows from (2.25) in view of the 
diffeomorphism invariance of the interaction term in the 
action. 



III. BIANCHI SPACETIMES 

In what follows, we shall assume both metrics to be 
homogeneous but anisotropic, that is, invariant under a 
three-parameter translation group G3 acting on the 3- 
space. The group is generated by three vector fields 
satisfying commutation relations 



(3.1) 



with constant structure coefficients C"^^j. Such groups 
have been all classified by Bianchi. The Jacobi identities 
imply that structure coefficients can be parameterized as 



(3.2) 



with n"'' — diag[n(-^), n^^), n^^)]. The different choices of 
the four parameters n^'^\ n^^\ a correspond to the 

nine Bianchi types. 

Denoting uj'^ the 1-forms dual to e^, the two metrics 
can be parameterized as 

-ait)^dt^ + habit) uj" (g>uj'', 

-Aitf dt^ + Habit) to" (E> Lo'' . 



dsi 



ds) 



(3.3) 



In this paper, we discuss only the class A Bianchi mod- 
els, for which the parameter o in (3.2) vanishes. These 
include types I, II, VIq, VIIq, VIII, IX. For these types 
tensors hab and "Hab can be chosen to be diagonal, and 
this guarantees that Goa — and Q^a — 0, so that there 
are no energy fluxes. To achieve a similar no-flux condi- 
tion for the type B (tilted) Bianchi classes with a 7^ 0, 
one has to either let hab and "Hab be non-diagonal, or to 
impose constraints on the otherwise independent values 
of their components. 

Let hab, Hab be diagonal matrices, 

hab = diag[Q;f,a|,a3]: 
■Hab = diag[^^ ^|,^|], 



and introduce 

0"6-2diag 

Q°b = 2diag 



ai a2 as 
ai ' as ' as 

Ai_ M A3 
Ai' A3' A3 



where the dot denotes the time derivative. The non- 
zero projections i?^ and TZ^g of the Ricci tensors on th« 
tetrad base = idt,u!°') and Ea = idt,ea) read 

[Q] , [Q^ 



pa 

n h 



\/hQ\ j (3) 



7^"o = + i4i, ^ — + n%. is a) 
° 2 4 ' " 2VnA 



Here the bracketed quantities are calculated according to 
(2.6), while the three-dimensional Ricci tensors 

(3) 1 

h R\ = 2iN^)\ - [N] N\ - i[N^] - -[NY)5\ , 



nn% = 2iN-')% - - iW'] - lwr)6% , (3.5) 



where 

N% - n'^'hcb, Af% = n-'^Ucb. 
The Ricci scalars are 



(3.6) 



(3) 1 

2h 



(3) 

R ^ ([iV]2 - 2[N^]) , n 



1 

m 



'[NY^2[N^]) 



(3.7) 

The non-trivial tetrad projections of the tensor 
1^^ = V¥^TZ are 



7 B = diag 



A Ai A2 A3 

a ' ai ' a2 ' as 



(3.8) 



Using this in Eqs.(2.4),(2.5),(2.20),(2.21) reveals that 
non-trivial tetrad components of the energy-momentum 
tensors t'''^'''^ and T^'^^^u are also diagonal. For example, 
the non-trivial tetrad projections of defined in (2.21) 
are 



T ^(no sum) = 

\a h + Xb + b3 AbAc 

B^A B.,C^A-B<C 

which defines the non-trivial components 



64A0A1A2A3, 



1 



T'^''^(no sum) = [r^(no sum) - "^j , 



T^'''''^(no sum) = — 



1 A^M 
a\/h 



T^(no sum) . 



with 



=hn + hi^\B + ■■■ + hi A0A1A2A3 . 
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This implies that 

j '^y^d*x = J{aUg+AUf)d'^x (3.9) 

where 

Ug = Vh{bn + h{Xi+X2 + X3) (3.10) 
+ 62(AiA2 + A2A3 + A3A1) + 63A1A2A3} (3.11) 

and is obtained from this by replacing h % and 
6fe — )■ hk+i- 

Finally, the matter energy-momentum tensors (2.23), 

assuming the fluid four-velocities to be tangent to the 



timelines and pg, Pg, pf, Pf to depend only on time, are 



T^"'^^^ =dmg[- Pg,Pg,Pg,Pg], 

r['"ll=diaghp/,P/,P^,P/]. 



(3.12) 



Inserting the above expressions (3.4)-(3.12) to the Ein- 
stein equations (2. 15), (2. 16) gives a system of non-linear 
ordinary differential equations for the field amplitudes 
a. A, aa, Aa, Pg, Pg, Pf, Pf which depend only on time. 

It is instructive to rederive these equations in a differ- 
ent way, by first inserting the metrics (3.3) to the action 
and then varying with respect to the field amplitudes. 



We adopt the following parametrization: 
The action (2.3) then assumes the form 



inn?S = — r / ae 
2kI 



1 



I fi / . . \ (3) . 



-tf + Pl + + R)d 



2.) 



Ae 



,3W 



A^ 



s>2 



(3) 

n 



\ f {aUg + AUf}d'x + 44'"' \9, Pa-Ps] + ^4"' [/' Pf^^S^ 



where 



Ug = 6oe'" + 636^^ + (^e-'iB+-M + 2e«+-''+ cosh[^(6_ - ;3_)]) 

+ 626^^+" (e2(«+-''+) + 26-^"+-^ cosh[^/3(S_ - /3_)]) , 
while Uf is obtained from this by replacing — !■ bk+i, whereas 

R = 2n(l)n(3)e-2f^e-2(/3+-^/3-) _ i g-^f^ |„(l)g2(^+ + y3/3_) _ „(2)g2(/J+-y30_) ^ ^(3)^-4/3+ 

(3) 

and TZ is obtained from this by replacing Q ^ W and f3± — >■ i3±. 
Varying the action with respect to a,fl,l3± yields the equations 



„3n 



n 



„3n 



n 



a 



„30 /5+ 



cos 77 



„3r2 



/3- 



(3) 

2 cos^ ?7 e^^Ug - e'^^R+ 2e^^Pg 

(3) 



9C/ 



3aUg) - 2ae^^R +'iae^^{pg ~ Pg) 



„30 P± 



^ ^ '2cos2r;[/-ae3"S I , 



i2a/3d 



„3W 



„3W 



,(3) 



while varying with respect to A, W, B± one obtains 

2 

2sin2r;e3^Z^/- e6^7^-H2e6W^/ 

\ (3) 

■iAUf \-2Ae^^n+ iA {pf - Pf) 



A 



„3W 



3W B± \ _ 



• 2 fdU 

sm 77 - — 



A 



1 d 



2sm^ rjU -Ae^^n 



(3) 



(3.13) 



d^x 



(3.14) 



(3.15) 



(3.16) 



(3.17) 
(3.18) 
(3.19) 

(3.20) 
(3.21) 
(3.22) 
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Here U = aUg + AUj. The matter densities and pressures satisfy the conservation conditions 

Pg+ih{pg + Pg)=0, f) f + 3W (f f + Pf) ^ . (3.23) 

DifFerentiating the first order constraint (3.17) and using (3. 18), (3. 19) to ehminate the second derivatives gives 



d 



d 



(3.24) 



which is nothing but the condition of conservation of T ^. DifFerentiating the first order constraint (3.20) and using 
(3. 21), (3. 22) to ehminate the second derivatives reproduces the same condition again, because T^'^^t is automatically 
conserved as soon as t''''''^ is conserved. 

Since the structure of the two line elements in (3.3) is invariant under time reparametrizations, a gauge condition 
can be imposed. For example, one can fix the gauge by requiring that a ~ 1. 

Equations (3.17)-(3.22) are equivalent to those obtained in the component approach by inserting (3.4)-(3.12) to 
(2. 15), (2. 16). The equivalence is seen in view of the relations 



(3) (3) 



(3) 



2R' 



(3) 

1 dR 



2 5/3+' 



(3) 

idR 



(3) (3) 



R\ = 



and 



2^m rj,[-i]a 





.[7]1 



(3) 

1 dR 

271 



(3.25) 



dU 



d/3_ ' 



(3.26) 



-[7]m 



(3) (3) M, 

as well as those obtained from these by replacing R% 11%, T'^^t V""t, Ug Uf, a A, n -i' W, P± -i' B±. 

Assuming the matter to consist of several components labeled by i with pressure being proportional to the energy 
density for each component (for example radiation plus a non-relativistic component), the matter densities and 

pressures determined by (3.23) are (with constant p^g^ and lUg*'') 



^WpWe-3(i+™«)n_ 



(3.27) 



Similar expressions for pf, Pg are obtained by replacing the index g ^ f and Q, — > W. 



IV. BIANCHI TYPE I SOLUTIONS 



A. Recovering General Relativity 



Let us first discuss the simplest Bianchi type 1 model, 
in which case the isometry group G3 acting on the 3-space 
is abelian. The two metrics are 

a=l,2,3 

ds}=~A{tfdt^+ J2 -^IWid^'f^ (4-1) 

a=l,2,3 

with tta, Aa parameterized according to (3.13). Since the 
spatial curvatures vanish, 

(3) (3) 

R% = TZ% = 0, 

the basic equations simplify, which allows us to find some 
exact solutions. 



If we assume the two metrics to be proportional, 

= C^g^... (4.2) 

then we obtain the GR solutions. Indeed, one has in this 
case 7^^^ — CS^^ and so 

r^, = ih + 362 C + 363 + bi C^)CS'', , (4.3) 

which gives the energy-momentum tensors 

«2rMt = -A,(C)<5^,, 

"1'^ --A/(C')<5^,, (4.4) 



with 



Ag{C) = cos^ 77 (60 + 36i C + 362 + h C^) 



sm rj 



(61 + + ibsC^ + b^C^ ) . (4.5) 
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Since the energy-momentum tensors should be conserved, 
it follows that C is a constant. As a result, we find two 
sets of Einstein equations: 



(4.6) 
(4.7) 



Since one has G"^ = G'^./C'^, it follows that A/ = Ag/C^ 
which gives an algebraic equation for C, 

cos^ (feo + 35i C + 862 + b3C^) 



sm 1] 
C 



(h + 362C + + 64C^ ) . (4.i 



It follows also that the matter sources should be such 
that 



As a result, the independent equations are (4.6), which 
are the same as in GR. In the present cosmological model, 
choosing the gauge where a = 1, these equations are 
obtained from (3. 17), (3. 18), (3.19), 



(4.10) 



where u± are integration constants. Solutions of these 
equations are reviewed in the Appendix. The solution 
for the f-metric is 



± • 



(4.11) 



As shown in the Appendix, if Ag is positive then the so- 
lutions approach the de Sitter metric exponentially fast, 
such that at late times one has 

17 = i/t + 0(e-3«*), ^± =^±(oo) + 0(e-3^^t), (4.12) 



with H = ^ Kgli. It is worth noting that, since we 
are in the Bianchi type I, the asymptotic values of the 
anisotropics, /3±(oo), can be set to zero via rescaling the 
spatial coordinates, which would not be true for other 
Bianchi types. 

If we choose the coupling constants according to 
Eq. (2.13), then Eq.(4.8) factorizes as 

(C7-1)P3(C^) = 0, (4.13) 

where P^(C^ is a cubic polynomial, 

Pz{.C) = (C3 + Ci)C^ + (3 - 5c3 + {x- 2)c4)Cf4.14) 
+((4 - 3x)c3 + (1 - 2x)c4 - 6)C + (3c3 + a ^ l)x, 

with X — tan^ ry, while 

Ag = cos^ 77(1 - C) (4.15) 
X ((c3 -I- C4)C^ + (3 - 5c3 - 2c4)C + Acs + C4 - 6). 



Depending on values of 03,04,77, the equation (4.13) can 
have up to four real roots. For example, for C3 = 1, 
C4 = 0.3, 77 = 1 the four roots are 

C = {-2.245; 1; 0.068; 3.616}, (4.16) 

and the corresponding 

Ag = {10.126; 0; -0.509; -4.505}. (4.17) 

As a result, we have four different solutions with four 
different values of the cosmological constant, which can 
be positive, negative, or zero. Choosing C = —2.245 
gives the de Sitter solution with the Hubble rate 

H = JKg/2, = 1.837. (4.18) 



(4.9) Below we shall find this value again for more complex 



solutions. If we choose C — 1, then Ag = Aj = and the 
two metrics and matter sources are identical: 



9tJ.v — ffiu , Pg ^ Pf J Pg ^ Pj 



(4.19) 



In vacuum, pg — Pg =0, the solution is either Minkowski 
metric if (t± = or the Kasner metric for non-zero a±. 



B. Solutions with equal anisotropies 

If the two metrics are different and the matter sources 
are not adjusted to be the same, then the simplest so- 
lutions of equations (3.17)-(3.22) are of the FLRW type 
[14], in which case anisotropies vanish. 



Pi 



B± = 0. 



(4.20) 



It turns out that one can obtain also more general solu- 
tions for which the anisotropies do not vanish but are the 
same in both sectors. 



3± = 6± ^ 0. 



(4.21) 



It is instructive to describe at the same time both 
isotropic solutions and solutions with equal anisotropies. 

The key point is that configurations with equal 
anisotropies extremize the potentials Ug and Uf so that 
their derivatives in (3. 19), (3. 22) vanish. 



df5± 

dUa 



.0 ^ 

/3±=B± ' 9/3± 



P±=B± 



0, 



dUf 



/5±=B± 



= 0, 



0. 



(4.22) 



As a result, Eqs.(3.19),(3.22) can be integrated to give 

„3n P± 



^3W ^± _ 5- 



a A 
with constant a±, S±. A particular choice 

a± =0, S±= 0, 



(4.23) 



(4.24) 
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corresponds to isotropic FLRW cosmologies, since 
Eq.(4.23) requires in this case that /3± = B± are con- 
stants, which can be set to zero by rescaling the spatial 
coordinates. 

If (7±,S± do not vanish then, since /3± = B±, taking 
the ratio of the two expressions in (4.23) gives the relation 



(4.25) 



As we shall see, the constant here should be positive, so 
that it is denoted by C^. 

In view of (4.22), the constraint (3.24) reduces to 



(4.26) 



which can be transformed to 



a I 



(4.27) 

Depending on which factor here vanishes, there are two 
solution branches, which we shall call generic and special. 



1. Generic solutions 

Let us first consider the case where the first factor in 
(4.27) vanishes, 



(4.28) 



Denoting 



one has 



w-n 



(4.29) 



The remaining equations to be solved are the constraints 
(3. 17), (3. 20), which reduce to 



-,2 2 i^2-6n I ^giO + Pg 



W =AU S'e 



(4.30) 
(4.31) 



Here a"^ = + cr^ and 5^ ^ S\ + S"^, while A^, A/ are 
obtained by replacing in Eq.(4.5) C — >■ ^. 

If we multiply (4.31) by {^a/AY and then subtract 
from (4.30) then, in view of (4.29), the result will be 



eKf{o-Kg{o^pg-ePf+£ 



(4.32) 



with 



(i) Isotropic solutions.- If ct^ = 5^ = then £ = 0. Since 
the matter densities determined by (3.27) are Pg{^) and 
Pf{W) with W being function of il and ^, Eq. (4.32) 
gives an algebraic relation between ^ and f2, which can 
be resolved to determine £,{ft). Inserting this function 
into (4.30) gives the equation for fi. 



,2 _ Agim) 



-Pgi^) _ i:.2 



(4.34) 



where we have set a — 1. Depending on the choice of the 
solution the form of the potential Ecs{ft) is differ- 

ent, leading either to self-acceleration or to other types 
of behavior [14]. If one requires the graviton mass contri- 





Eeff 






\^^^^^^^ 4=exp(W-n) 







£2 



FIG. 1. The solution ^{fl) of the algebraic equation (4.32) 
for C3 = 1, C4 = 0.3, 77 = 1, pg = 0.25 x e"*" + 0.25 x e"^", 
Pf =0, and the corresponding potential i5off(f2) in (4.34). 



bution to the total energy density, A^, to be dominant at 
late times (accelerating phase) but small at early times 
(matter dominated phase), then one has to have [14] 



as well as 



r,w in 



as 



&i < 0. 



0, 



(4.35) 



(4.36) 



The function ^(51) is then always negative'^. 

An example of a self-accelerating solution is shown in 
Fig.l. For il —i' 00 one has 



-2.245, E^si^) H ^ 1.837, (4.37) 



where H is the Hubble expansion rate. In Section V we 
shall study more general solutions with exactly the same 
late-time behavior. 

(a) Anisotropic solutions.- If = a^C^ 7^ 0, then com- 
bining (4.28) with (4.25) gives 



(4.38) 



1 - 



(4.33) 



^ f can be negative, since or need not to be positive definite, 
because only their squares appear in the metric coeflicients. 
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from which, denoting the integration constant by v, we 
find 



C 



(4.39) 



If V does not vanish, then 

£ = ■iva'^e-^'^ (2 - ve^'^) . 

Inserting this to (4.32) gives an algebraic relation be- 
tween ^ and f2, which, in view of (4.39), becomes a rela- 
tion between and C. As a result, 17 should be constant 
and no dynamical solutions exist in this case. 

If 1/ = then i = G and Eq.(4.25) gives Ca%Q 
that the two metrics are proportional. 



(4.40) 



and £ = 0. This case was discussed above in Section 
IV A, where we saw that the matter sources should be 
fine-tuned, Pg = C^Pf- 

2. Special solutions 

Let us now consider the case where the second factor 
in (4.27) vanishes. 



he + 262^ + 6i - 0, 



(4.41) 



implying that ^ is constant. Since it should be real, one 
has to have h"^ — > 0. The equations to be solved 
are again (4.30) and (4.31), where now W — Cl. One has 
either cr^ = 5^ = in the isotropic case or. in view of 
(4.25), 



(4.42) 



in the anisotropic case. In both cases, taking the differ- 
ence of (4.30) and (4.31) gives 



A 



AfiO+Pf 



a . 



(4.43) 



In the isotropic case this completes the procedure - so- 
lution for g^j^iy is the same as in General Relativity with 
the cosmological constant Ag(^) and matter pg. The sec- 
ond metric is obtained using the values of = ^e^ and 
of A from (4.41), (4.43). The solution exists if only the 
argument of the square root in (4.43) is positive, which 
condition turns out to be rather restrictive. For example, 
if the parameters bk are chosen according (2.13), then one 
finds 



COS^ T] 



sin^ rj 



(4.44) 



so that if Ag > then Af < 0, therefore only solutions 
with pf > ~Af are allowed. 



For anisotropic solutions Eqs.(4.42) and (4.43) give 



Agio + Pa 



(4.45) 



The argument of the square root here should be positive 
and, in addition, constant. Since, by our assumption, 
Pg, Pf do not contain constant contributions, they should 
be proportional, pf — {Af / Ag)pg. This excludes self- 
accelerating solutions with Ag > 0, since if Ag and pg 
are positive, then both Af and pf are negative and the 
argument of the square root is negative. 



V. MORE GENERAL BIANCHI TYPES 

For the general Bianchi class A models we choose 



a=l,2,3 



ds} ^ ~ A{tfdt^ + J2 Alit) uj" (E) iv" , (5.1) 



a=1.2,3 



where the vectors dual to uj"- do not commute, so that 
the spatial curvature do not vanish. Few solutions then 
can be found in a closed form. 



1. Recovering General Relativity 

Assuming again the two metrics and their sources to be 
proportional, /^^ ~ C^g^^,, we obtain the GR equations. 



f]2 



Pl+Pl-\.R + '^, (5.2) 
6 3 

(3) 

e^" dR 



(5.3) 



(3) 



with R given by (3.16) and Ag,C defined by (4.5), (4.8). 
One can choose C = 1, in which case Ag — 0. 

One cannot integrate the equations for /3± analytically, 
unless one assumes that /3± = 0, which is however consis- 
tent with (5.3) only when n^^^ = n^"^^ = nP'^ = k, which 
corresponds to the Bianchi I,IX for fc = 0, 1, respectively. 
Eq.(5.2) then becomes 



-20 



A„ 



(5.4) 



which describes the spatially fiat {k = 0) or spatially 
closed {k — 1) FLRW universe. 



2. The isotropic case 

If the the anisotropics vanish but the metrics and their 
sources are not proportional, then instead of (5.4) one 
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obtains 



(5.5) 



where i?cff(f7) is the same as in (4.34). It is instructive 
to rewrite this equation as 

a? + V{a) ^ -k (5.6) 

with a = 2e" and V = -s?El^- Since a? > 0, the 
solution is restricted to the regions where 

V{a.)<-k. (5.7) 

To construct ^(a) one should resolve the algebraic equa- 




n=ln(a/2) 

FIG. 2. Potential V^(a) in (5.6) for isotropic solutions for 
different choices of the matter density: pi = 0.1 x e~*^ and 
P2 = 0.25 X e""" + 0.25 x e"^" and pf = 0. In the Bianchi IX 
case the motion is restricted to the region V < —1, so that 
for pg = pi one has a recoUapsing solution corresponding to 
the region on the left from the barrier, and a bounce solution 
corresponding to the region on the right. 

tion (4.32) (with £ = 0) to obtain ^(il) for chosen pg, pf, 
and then insert the result to (4.34) to obtain Ecs{^). In 
Fig. 2 the result is shown for two different choices of pg 
and ioT Pf — 0. One can see that, since V < 0, for k — 
(Bianchi I) the motion covers the whole region a > 0, 
including the initial singularity at a = 0. The universe 
therefore expands forever starting from zero size. 

The same is true for k — 1 (Bianchi IX), but if only 
Pg is large enough (pg = p2), so that the total amount 
of matter in the universe is sufficient. If pg is small 
{pg — pi) then the total 'energy' is not enough to pass 
over the barrier, and the motion is divided into two re- 
gions. The region on the left from the barrier corresponds 
to universes which start at the singularity, then expand 
up to a maximal size, reflect from the barrier, and recol- 
lapse. The region on the right from the barrier describe 
universes that shrink from inflnity down to a minimal 
non-zero size, then bounce from the barrier and expand 
again to infinity^. As we shall see later, these features 
apply also when the anisotropics are taken into account. 



A. Generic solutions — integration procedure 

As we have seen, exact solutions can be found either 
in the isotropic case of the Bianchi I and IX, or when the 
anisotropics are equal, uniquely for the Bianchi I. In all 
other cases we are bound to use the numerical analysis. 
Our proceure will be described below, while its outcome 
can be summarized already now: it turns out that the 
solutions with equal anisotropies in the two sectors play 
the role of asymptotic states to which the spacetime ap- 
proaches at late times. 

In order to integrate the equations numerically, we flrst 
of all convert the equations to the form of a dynamical 
system. Let us introduce the variables 



y3 

ye 
y9 



„30 



a 

o3W 



2/1 = 

y4 -- 

2/7 



2/2 



(5.8) 



o3f2 



„3n 



2/8 



■I3-, 



o3W 



o3W 



A 



2/10 



A 



B+, 2/11 



A 



The second order field equations (3.18), (3.19), (3.21), 
(3.22) can be represented in the first order form 



Vn = FN{a,A,yM) , 



(5.9) 



where iV, M = 0, . . . , 11 and F/v(q!, A, yn) are defined in 
the Appendix. 

The first order equations (3.17), (3.20), (3.24) are con- 
straints 



Ci(2/jv) 
with 

Ci = -vl 

C2 



0, C2(yAr)=0, C-i{yN)^0, 



(5.10) 



2/? 



vl 



1 (3) 

2^vlUg--ylR 



72/0 Pg> 



yl 



1 (3) 1 

Vw + 2/n + 2s2y3z^^ ^_yln + -ylpf, 



2/3 2/62/0 



2/0 2/92/3 



d 
dyo 
d 



d 

2/72/1 

oyi 
d 

yioyijr- 

02/4 



VSyii 



2/5 



dys 



(3) (3) 

where Ug,Uf, R, TZ, pg, pj are defined in the Appendix. 

In order to implement the constraints and also deter- 
mine the lapses a and A, we remember that the third con- 
straint in (5.10) is obtained as the condition of propaga- 
tion of the first two. Specifically, calculating the deriva- 
tives 



Ci 



N=0 " 



V F 

^dy^ " 

N=0 " 



(5.11) 



* The similar bounce is found in de Sitter solution for the closed 
universe, i.e., a oc cosh(Ht). 



gives expressions proportional to C3. Therefore, if the 
third constraint is fulfilled, C3 = 0, then Ci = C2 = 0, 
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so that it is sufficient to require Ci = C2 = only at the 
initial moment of time. 

Now, what guarantees that the third constraint prop- 
agates itself ? In fact, calculating its time derivative 

c,.i:gF„, (5,12) 

and setting the result to zero we discover a non-trivial 
relation of the form 

aXa+AXA^O, (5.13) 

where Xa and X^ are rather complicated functions of 
UN (we do not show them explicitly). This is the condi- 
tion of propagation of the third constraint, which in turn 
guarantees that the first two constraints propagate. This 
condition determines the lapse A, 

^=-|^a. (5.14) 

Our procedure is then as follows: we integrate the 12 
equations (5.9) with A determined by (5.14), where we 
can use the gauge condition 

a = 1. (5.15) 

The three constraints are imposed at the initial time mo- 
ment, and the above consideration guarantee that they 
are fulfilled for all times. 

To determine the initial data, we choose 9 out of 12 
amplitudes y^r and give them some initial values. It is 
convenient to choose 

yo, yi, y2, y4, yd, y?, ys, yio, yn (5.16) 

which determine O, f3±, B±, $±, B± at t = 0. Having 
chosen these 9 values, the values of y^^y^, yg determining 
W, (l, W are obtained by numerically solving the three 
constraint equations (5.10). Having done these, all ini- 
tial data are determined, and there remains just to inte- 
grate forward the 12 equations (5.9) using the numerical 
routines of [19]. 

To check our procedure, we specialize to the Bianchi I 
type by setting in the equations n^^^^ = rt^^^ = n^'^^ = 
and also switch off the matter, pg = pf = 0. We then 
first consider the case where all anisotropics and their 
time derivatives initially vanish, 

2/1 = 2/2 = y4 = = 1 , y? = 2/8 = yio = yu = 0. 

The numerical integration then reproduces the de Sitter 
solution, for which these conditions are preserved for all 
times. Next, if we initially set 

yi = y4 7^ 1, y2 = ys 7^ 1, y? = ys = yio = yn = 0, 

then these conditions are also preserved for all times. 
This corresponds to the solution with a± = S± = but 



with constant values of /?+ = B+ and /3_ = This is 
again the de Sitter solution, since constant anisotropics 
can be removed by rescaling the spatial coordinates. 
As the next step, we initially choose 

yi=y4 7^l, y2 = y5 7^i, y7 = o'+, ys = <T^, 
yio = C<7+, yii=Ca_, C^yl/{Ayl). (5.17) 

The numerical solution then again preserves these con- 
ditions at all times, which corresponds to the analytic 
solutions with equal anisotropics. 

B. Numerical results 

After consistency checks, we turn to the generic case. 
We switch on the matter by setting in the equations pg ^ 
and Pf 0. We choose an initial value for yo = and 

some values for yi — e^+, ?/2 — e^^- , y^ — e'^+ , 2/4 = 
g>/3B_ assuming that f3± and B± are not too large. In 
addition, we choose small enough values for j/y, y^, yiQ, yn 
determining /3±, B±, so that the initial configuration is 
not too anisotropic. 

We consider all Bianchi class A types. The values of 
the parameters n^"^ are given in Table I. 





I 


II 


VIo 


VIIo 


VIII 


IX 







1 


1 


1 


1 


1 










-1 


1 


1 


1 
















-1 


1 



TABLE I. Values of n^"' for Bianchi class A types 

We find that if the parameters are chosen such that the 
universe enters the self-accelerating regime, then generic 
solutions rapidly approach the state with equal and con- 
stant anisotropics, /?+ = B+ and /3_ = B-. For Bianchi I 
type one can scale away these constant values by rescaling 
the spatial coordinates, but this cannot be done for other 
Bianchi types. As a result, the final state of the universe 
is anisotropic, even though the metric amplitudes Cl, W, 
A approach constant values, as for the isotropic universe. 

Some examples of our solutions are presented in Figs. 3- 
9 for the parameter values C3 = 1, C4 = 0.3, ij = 1, which 
are the same as for the FLRW self-accelerating solution 
in Fig.l. We assume the f-sector to be empty, while the 
g-sector to contain a radiation and a non-relativistic mat- 
ter described by pg = 0.25 x e~^^ -1-0.25 x e^'^^ , which is 
the same function as in Figs. 1,2. There is nothing special 
about this choice, since for other parameter values solu- 
tions are qualitatively the same. (We notice that ~ 1 
for f7 = 0, therefore, if m ^ 10^'^'^eV, the dimensionful 
energy density in (2.23) is ~ vr? j ~ 10~^°(eV)^, which 
is about the present density of the universe). The uni- 
verse will accelerate if C3 , C4 are chosen such that 61 < 
(see Eq.(4.36)) while pg is large enough in order that the 
system could travel over the potential barrier as shown 
in Fig. 2. It is also worth noting that the time scale in 
Figs.3-9 is the Hubble time, l/H. 
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FIG. 3. The Hubble parameter Q. for solution for all Bianchi 
types of class A. 
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FIG. 4. Anisotropy amplitudes profiles for Bianchi IX. For 
other Bianchi types the picture is similar. 



The initial value of the scale factor is 17 = and the 
initial anisotropics are /3+ = 0.02, /3_ — —0.013, i3+ = 
—0.02, B- = 0.015, while the initial values of B± are 
set to zero. Again, changing these values does not change 
the qualitative behavior of the solutions. The values of 
yV, (l, yV are determined by the constraints (5.10). 




0.5 1 1.5 



Ht 

FIG. 5. The ratio of the two scale factors ^ = e^"*^ and the 
lapse amplitude A for all Bianchi types. 



As seen in Fig. 3, for all Bianchi types the universe 
expansion rate rapidly approaches the constant value, 
Cl ^ H = 1.837 (the same value as in Eq.(4.18)), so that 
the universe approaches the de Sitter phase. If we return 
for a moment to dimensionful quantities, the Hubble rate 
becomes H — Hm, so that the cosmological constant is 
indeed due to the graviton mass. 
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FIG. 6. The relative magnitude of shears E = 0l+l3iy^'^/Q 
for all Bianchi types. 

The anisotropy amplitudes f3±, B± approach con- 
stant values within an approximately twice Hubble time, 
2H^^, as is seen in Fig. 4, where the solutions for the 
Bianchi IX are shown. The ratio of the two scale factors 

= e^/e^ is shown in Fig. 5. It is always negative for all 
solutions and rapidly approaches the value C — —2.245, 
which is the same as in Eq.(4.16). The lapse function 
A, also shown in Fig. 6, approaches the same value, so 
that the two metrics become proportional at late times, 
ffiv = C^gtj.u, as for the solutions in section IV A. 

Fig. 6 shows the relative magnitude of shears 

^ ' 

which measures the relative contribution of the 
anisotropics to the total expansion rate. The maximal 
value of E depends on the initial anisotropy values. How- 
ever, setting the latter to zero we obtain practically the 
same curves for the Bianchi types II,VIo,VIIo,VIII, since 
anisotropics cannot be zero in these cases, so that they 
are driven to non-zero values by the cosmic expansion. 
The shears approach zero exponentially fast in Hubble 
units. However, if our universe entered the acceleration 
phase only recently, then one can expect that only few 
e-folding times have elapsed since then, so that the shear 
energy is not necessarily small at present. 

The shears /3±, B± multiplied by e^^^^ oscillate with 
constant amplitudes, as shown in Fig. 7, so that $± ~ 
B± ^ e"'^^/^. To understand this, we remember that our 
configurations approach the solution with proportional 
metrics, f^^^ = C^g^,j, described in section IV A, so that 
the deviations from this solutions are small at late times. 
We therefore linearize the field equations with respect to 
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Ht 

FIG. 7. The rescaled shears e^^^'^P± (smaller amplitudes) 
and e^^^'^B± (larger amplitudes) for the Bianchi IX solution, 
for other Bianchi types the picture is similar. 

the deviations, and solving the linearized equations we 
obtain 

/3± ~ 6± ~ e-^"^'^ cos{Hujt) (5.18) 

with 

2 /, f,^, ^2, .fCcos^-q sin^?7\ 9 
- = {hi + 2Cb, + C%) [—^ + 

(5.19) 

For our parameter values this gives uj — 1.183 and the os- 
cillation period in Hubble time units T — 2tt/uj — 5.309, 
in perfect agreement with what is shown in Fig. 7. 

Therefore, the shear contribution to the Hubble rate is 

I3l+I3^_^e-^^^ r^l/a.^, (5.20) 

which falls off similarly to the energy density of a non- 
relativistic matter. In GR (see the Appendix) the shear 
fallofF is much faster, + ^ 1/a^, corresponding to 
a 'stiff matter'. Since in the bigravity the shear contri- 
bution to the total energy density increases slowly, this 
could have observational effects. 



1. Near singularity behavior 

It is interesting to see how the solutions continue to the 
past. Their parameters are chosen to avoid the bounce 
behavior, so that when continued to the negative t re- 
gion, they should hit a singularity at some point. The 
numerical simulations confirm these expectations and re- 
veal that for all Bianchi types under consideration there 
is a singularity for t < 0, where both and vanish. 
For the Bianchi I this is shown in Fig. 8. Interestingly, in 
the negative t region the solution shows a throat, and 
first expands before collapsing to zero. 

Let us consider in more detail the case of Bianchi IX. 
In General Relativity Bianchi IX solutions reveal chaotic 
features. When approaching singularity, the metric co- 
efficients aa which measure the proper distances along 
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FIG. 8. The scale factors near the singularity. 

the spatial axes (defined in (3.13)) show an infinite num- 
ber of oscillations, whose positions and amplitudes are 
ergodic [20]. Within an effective description, such a be- 
havior is explained by a two dimensional 'billiard' motion 
of a particle reflecting from rigid walls. 
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FIG. 9. The metric coefficients aa, Aa near the singularity. 

At first glance, nothing similar is seen in our case. 
Fig. 9 shows the Bianchi IX solution continued to the 
past, and the singularity corresponds to a point where 
both and vanish, but in its vicinity coefficients 
aa, Aa defined by (3.13) approach zero without oscilla- 
tions. However, zooming the picture one can see that os- 
cillations actually start, but not many of them are seen, 
since we cannot approach the singularity close enough 
due to numerical errors. 

The picture becomes more clear for solutions with- 
out matter. Setting pg = Pf — 0, the simplest salf- 
accelerating solution is pure de Sitter, but it is of course 
non-singular. It can be obtained from Eq.(5.6), which 
reduces to (remember that a = 2e^) 

s?=H'^s?-l a= 4 cosh(t-io), (5.21) 

H 

where H = 1.837 is the same as for all other solutions 
under considerations. Qualitatively, the singularity is 
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FIG. 10. Self-accelerating for t > Bianchi IX solutions 
continued to the t < region. Choosing at t = /9± ~ 
B± ~ 10~^, the solution is of regular bounce type, but if 
j3± ~ B± ~ 10"^ then the solution develops a singularity 
where vanishes. 

avoided because a cannot be too small, since otherwise 
the right hand side of the equation would be negative. 
Therefore, a bounces back when it achieves the minimal 
value 1/H. 

However, for non-zero anisotropy the equation in (5.21) 
receives in the right hand side additional terms propor- 
tional to -|- , and these can keep the whole expres- 
sion on the right positive even if a — > 0. Therefore, so- 
lutions with high enough anisotropy can approach sin- 
gularity. To verify this, we choose at i = two sets of 
initial values for the anisotropics, f3± ^ B± ^ 10^^ and 
P± ~ S± ~ 10~^. When continued to the future, we see 
self-acceleration in both cases, but when continued to the 
past, we obtain the bounce behavior in the first case and 
a curvature singularity in the second case (see Fig. 10). 
Therefore, the empty de Sitter spacetime becomes singu- 
lar when too much anisotropy is added. 

It happens that for the empty Bianchi IX solution we 
can approach the singularity much closer numerically, 
and in Fig. 11 we show In(aa) against f2 cx ln(t). This 
time one can clearly see the typical features of the bil- 
liard motion characterized by a sequence of Kasner-like 
periods. During each period one has aa on , where 
Pa fulfill Eq.(A.l.lO). During the next period, the val- 
ues Pa change such that one of then remains positive, 
the corresponding amplitude continues to decrease to- 
ward singularity {a^ in the Fig. 10), while two other pa's 
change sign, such that the increasing amplitude becomes 
decreasing and vise versa. Of course, within our numer- 
ical approach we can capture only the beginning of the 
infinite sequence of Kasner-like steps, and so it would be 
interesting to develop an affective analytical description. 
If the mater is present, provided that w < 1^, then its 
contribution to the equations is subleading as compared 




FIG. 11. The anisotropy parameters Oa for the empty Bianchi 
IX near the singularity plotted against fi. They show the typ- 
ical 'billiard' behavior characterized by a sequence of Kasner 
periods. For Aa the picture is similar. 

to that of the anisotropy terms, so that it cannot influ- 
ence the chaotic behavior. Therefore, since the empty 
Bianchi IX solutions are chaotic, so should be those with 
matter. 



VI. CONCLUDING REMARKS 

We have studied anisotropic cosmologies in the ghost- 
free bigravity assuming both metrics to be diagonal and 
of the same Bianchi type of class A. Including a source 
consisting of a radiation and a non-relativistic matter, we 
considered generic initial data describing (not too large) 
anisotropic deformations of a FLRW universe. We find 
that the universe evolves into a state in which it ex- 
pands with a constant Hubble rate proportional to the 
graviton mass, while the anisotropy parameters approach 
constant non-zero values. In the Bianchi I case constant 
anisotropics can be scaled away be redefining the spa- 
tial coordinates, but not for other Bianchi types. For 
example, for the Bianchi IX solutions the constant t spa- 
tial sections will be not round 3-spheres but squashed 
spheres. 

The conclusion is that generic self-accelerating cos- 
mologies in bigravity are anisotropic. The anisotropy 
contribution to the total expansion rate approaches zero 
exponentially fast in Hubble units, but if our universe en- 
tered the acceleration phase only recently, then one can 
expect that only few e-folding times have elapsed since 
then, so that the shear energy is not necessarily small. At 
late times the shear contribution to the total expansion 
rate decreases much slower than in GR, only as an in- 
verse cube of the size of the universe, which is the same 
falloff rate as for a non-relativistic matter. Therefore, 



^ The chaotic behavior in Bianchi IX disappears if a massless scalar 
field is included, since the efi^ective equation of state in this case 



is P = p, so that near singularity p oc e ° is able to compete 
with the effect of anisotropics [21], 
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the anisotropy effect could be visible, although compar- 
ing with observations goes beyond the scope of this paper 
(see [22] for a data-fitting for the FLRW cosmologies with 
massive gravitons). 

The fact that the anisotropy contribution shows the 
same falloff rate as a cold dark matter suggests that the 
latter could in fact be the effect of the anisotropics. If 
this is true, then it follows that theories with massive 
gravitons can explain both the dark energy - the cosmo- 
logical term mimicked by the graviton mass, and also the 
cold dark matter mimicked by the anisotropics. At the 
same time, it is unclear if this interpretation can explain 
also the other properties of dark matter, as for example 
clustering around galaxies. 

It would be interesting to develop an analytic descrip- 
tion for the behavior of the Bianchi IX solutions near 
singularity. Yet another interesting open issue would be 
to see if the FLRW solutions with non-diagonal metrics 
studied in [8], [9], [10], ]11], ]12] could be generalized for 
non-zero anisotropics. 



The last of Eqs. (A. 1.2) can be integrated to give 

e3"/3± = a± , (A. 1.4) 

where a± are integration constants, from where 

P± = P±{to) + a± f e-^^'dt. (A.1.5) 



The equations then reduce to 

2 



(A.1.6) 



where and p is determined by the energy 

conservation condition. 



assuming the equation of state P — wp. 



(A.1.7) 



a. Vacuum solutions 
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Appendix A: Bianchi I spacetimes in GR 

In this Appendix, we discuss the Bianchi I spacetime 
in General Relativity for some matter types. Assuming 
the metric form 

ds^ = -df + J2 oia{tf{dx''f (A.1.1) 



a=l,2,3 



where 



the Einstein equations reduce to 

(e'^h) ^\e^'\p~P), 
(e3"/3±)' = 0, (A.1.2) 

where p and P are the energy density and pressure, which 
fulfill the energy conservation condition. 



p + itl[p + P)^Q. 



(A.1.3) 



Setting in the above formulas p = po = 0, one can 
integrate Eqs. (A. 1.5), (A. 1.6) to obtain 



Ua 'x.t^" , a — 1, 2, 3 



(A.l. 



where 



1 (T+ + \/3cr_ 1 cr+-\/3CT_, 

+ )' P2 = x(l + , 

6 (7 OCT 



P3 = ^(l-2^), 
o a 

so that 



(A.1.9) 



Pi+P2+P3=pI+pI+pI = 1- (A.1.10) 
This corresponds to the Kasner solution. 



b. Cosmological constant 

Let us choose in (A.1.7) w = —1 and p = po = 3H^. 
The solution of (A.1.6) is then 



Am 



(A.1.11) 



Inserting this to (A.1.5) determines If ctj- = then 
the solution is pure de Sitter, = H(t — to), with con- 
stant anisotropy parameters P± which can be set to zero 
via rescaling the spatial coordinates a;°. If ct± 7^ 0, then 
the solution approaches the de Sitter metric exponen- 
tially fast, 

n^H{t- to) + 0(e-3ff*), /3± = /3±(oo) + 0(e-3^*) . 

(A.l. 12) 

Calculating the integral in (A.1.5) gives explicitly 



A4 



a = 1,2, 3, (A.l. 13) 
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where pa are the same as in (A. 1.9) and 



2H 



(A.1.14) 



One has X± — > 1 as i — > oo, so that solutions approach 
de Sitter metric. On the other hand, taking the Hmit 
H and choosing the integration constant such 
that CT/(2ff)e3^*" = 1 one has X+ ^- 2 and X_ 3Ht 
in this hmit. and therefore aa oc t^" , so that the Kasner 
solution is recovered. 

One can similarly obtain solutions for a negative cos- 
mological constant, when w = —1 and p = po = —3H^. 
The solution of (A. 1.6) is then 



H 



sin(3iJ(i - to)) . 



(A.1.15) 



inserting which to (A. 1.5) gives 



3i/(i-<o)\'^V. 3Hit-to) 
an oc COS tan 



where Pa are again the same as in (A. 1.9). 



(A.1.16) 



which shows that the anisotropy contribution always be- 
comes small for large fi, provided that w < 1. This 
equation can be integrated in quadratures, the late time 
behavior of the solution is 



/3 = /3(oo) + 0(t^^). 



(A.1.18) 



One can also consider a more general matter consisting 
of several components with different equations of state. 
Pi = WiPi, in which case 

p = ^p,oe-3(i+».)n. (A.1.19) 



For example, if there is a cosmological constant plus a 
one-component perfect fluid, then Eq.(A.1.6) reduces to 



jj2^6n ^ Po^snn-^) ^ (A.1.20) 
3 



If w > —1 then for large Q the second term on the 
right becomes dominant and the universe approaches the 
de Sitter state. The conclusion is that in all cases the 
anisotropy effect soon becomes negligible when the uni- 
verse expands. 

Appendix B: Dynamical system description 



c. More general matter 
For a general equation of state Eq.(A.1.6) reduces to 
(e^f^fj)' =a^ + P^emi'^) ^ (A.1.17) 



The second order equations (3.18), (3.19), (3.21), 
(3.22) can be represented in the first order form 

i/N = FN{a,A,yM), (A.2.1) 
where the variables yN{t) are defined in (5.8), while 



t<o — a^, ti—a — o- , i^2 = v3( 



vl 



Vo 



Fa = 



2 I dU \ ,(3) 

cos '? ( 2/0^ + iaUgj + 3a y^ {pg - Pg) -2ayo R 



, F7 = 



12 



„ 2 9U 

-2 cos vyiT^ + yi 



(3) 

dR 
dyi 



12 



12 



(3) 

2 dU n dR 

-2cos riy2^ Vayoy2^ 

oy2 oy2 



(3) 

„ . 2 dU , ^ dR 

-2sm 772/4- VAy^yi^ 

dyi dyi 



Fa 



Fn 



sin' ^ ( 2/3^ + 3AUf] + 3Ayl (pf - Pf) -2Ayln 



73 
12 



(3) 

^ ■ 2 9U , dR 

-2sm 7jy5- hAy^ys^ 

dy5 dy5 



Here U = aUg + AUf with 

Ug = ylibo + bl{M + A2 + A3) + &2(AlA2 + A1A3 + A2A3) + 63A1A2A3} 

and Uf is obtained from Ug by replacing bk — )• 6fc+i, 



where the eigenvalues A^ are given by 



Ai = 



ysyiVb 
yoyiy2 



A2 - 



yoyiys' 



A3 = 



The 3-curvatures are 

Ml (A22) 1 f n(') y^f - n(^) yf + n^^) yl 

yoyl ■ ^ ■ ■ ' ylyl 2 yoyiz/i 
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. (3) ^ 

while TZ is obtained by replacing in this expression yo — > 
2/3j yi ~^ 2/4j ?/2 ~^ 2/5- Finally, the matter terms are 

i i 

(A.2.3) 



where Pg , Wg are constant parameters, while Pf,Pf are 
obtained by replacing in these expressions yg — > y^ and 
Pg-* — >■ py'', Wg*-* — >■ w^*'' (the number of matter compo- 
nents needs not to be the same in both sectors). 
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